Fingerprinting Chemical Markers in the Mediterranean Orange Blossom Honey: UHPLC-HRMS Metabolomics Study Integrating Melissopalynological Analysis, GC-MS and HPLC-PDA-ESI/MS

(1) Background: Citrus honey constitutes a unique monofloral honey characterized by a distinctive aroma and unique taste. The non-targeted chemical analysis can provide pivotal information on chemical markers that differentiate honey based on its geographical and botanical origin. (2) Methods: Within the PRIMA project “PLANT-B”, a metabolomics workflow was established to unveil potential chemical markers of orange blossom honey produced in case study areas of Egypt, Italy, and Greece. In some of these areas, aromatic medicinal plants were cultivated to enhance biodiversity and attract pollinators. The non-targeted chemical analysis and metabolomics were conducted using ultra-high-performance liquid chromatography high-resolution mass spectrometry (UHPLC-HRMS). (3) Results: Forty compounds were disclosed as potential chemical markers, enabling the differentiation of the three orange blossom honeys according to geographical origin. Italian honey showed a preponderance of flavonoids, while in Greek honey, terpenoids and iridoids were more abundant than flavonoids, except for hesperidin. In Egyptian honey, suberic acid and a fatty acid ester derivative emerged as chemical markers. New, for honey, furan derivatives were identified using GC-MS in Greek samples. (4) Conclusions: The application of UHPLC-HRMS metabolomics combined with an elaborate melissopalynological analysis managed to unveil several potential markers of Mediterranean citrus honey potentially associated with citrus crop varieties and the local indigenous flora.

In order to investigate the differences between the honey extracts, taking into account their geographical origin, a statistical analysis was performed based on the results from the UHPLC-HRMS analyses. Firstly, a heatmap was generated for both negative and positive modes to investigate the relative content of secondary metabolites and compounds among the different countries. As shown in Figure 3, citrus honey from Italy exhibited more differentially increased compounds in the comparisons (between countries).
In this context, PCA and OPLS-DA were calculated using R for visualization of any metabolic clustering of the different groups of samples. As presented in Figures 4 and 5, statistical analysis managed to separate honey samples when these were compared bilaterally (country level), in the negative mode for all countries ( Figure 4) and for Greece vs. Italy and Italy vs. Egypt, in the positive ion mode ( Figure 5). The statistical analysis revealed several compounds to be differentially increased among the three countries. These chemical markers are presented in Table 1 and emphasize the importance of both geographical and botanical origin, also considering other honeybee foraging plants apart  In order to investigate the differences between the honey extracts, taking into account their geographical origin, a statistical analysis was performed based on the results from the UHPLC-HRMS analyses. Firstly, a heatmap was generated for both negative and positive modes to investigate the relative content of secondary metabolites and compounds among the different countries. As shown in Figure 3, citrus honey from Italy exhibited more differentially increased compounds in the comparisons (between countries).
In this context, PCA and OPLS-DA were calculated using R for visualization of any metabolic clustering of the different groups of samples. As presented in Figures 4 and 5, statistical analysis managed to separate honey samples when these were compared bilaterally (country level), in the negative mode for all countries ( Figure 4) and for Greece vs. Italy and Italy vs. Egypt, in the positive ion mode ( Figure 5). The statistical analysis revealed several compounds to be differentially increased among the three countries. These chemical markers are presented in Table 1 and emphasize the importance of both geographical and botanical origin, also considering other honeybee foraging plants apart from installed AMPs. It is important to mention that the integration of the four non-citrus Egyptian honey samples in the metabolomics provided comparable findings when solely citrus honey was considered. Consequently, they were not excluded, and the totality of samples is presented. In the same context, one honey obtained from the Italian market (sample number 20) did not correspond to its label as citrus honey, verified both by HRMS and melissopalynological analysis.      As presented in Table 1, the honey samples from Italy compared to Greece and Egypt were differentially increased in glycerophospholipids, flavonoids, and some terpenic molecules. In the same context, in Egypt (when matched to Greece) fatty acyl glycosides of mono-and disaccharides and carboxylic acids were found to increase. Greece compared to Italy showed increased levels of selected terpenoids, iridoid glycosides, and fatty acids. These differences, presented in detail below, are potentially attributed to the different orange tree varieties used in the case studies of the project, especially between Greece and Italy, where the majority of citrus-based honey came from, and the geographical origin. A recent study showed the effect of geographical origin on the differences concerning the chemical composition among the same types of honey [33]. In another domain, recent work on the chemical composition of juices originating from several European pear cultivars has shown that such differences are cultivar-dependent [34].

Flavonoids
The use of flavonoids as chemical markers of honey's botanical origin has been reported by several research groups [35,36]. Similar patterns have been used to identify bee pollen's respective origin [37]. In this work, attention was given to discriminate citrus honey from three distinct Mediterranean countries. This task was challenging since the majority of samples were of citrus crop origin and the margins of differentiation among them, at first reading, seemed limited.
Hesperidin (flavanone) and its aglycone form, hesperetin, are characteristic flavonoids of citrus honey. Hence, insight into the variability of concentrations among the same type of honey is of the utmost importance. The latter becomes even more significant when attempts are made to discriminate honeys based on geographical criteria. It is noteworthy that hesperetin (tR = 12.71, with 302.0790 m/z Table 1) was increased in Italian honeys, while hesperidin (tR = 10.14 with 610.1903 m/z Table 1), a glycoside of hesperetin, was comparatively increased in Greek honeys ( Figure 6). Both molecules can be further investigated under a biosynthetic framework. Flavonoids' biosynthesis has been studied by many research groups, though for several of its members (i.e., hesperetin), it is still unknown. Hence, any observation of their formation and the conditions that govern them is significant. As presented in Table 1, the honey samples from Italy compared to Greece and Egypt were differentially increased in glycerophospholipids, flavonoids, and some terpenic molecules. In the same context, in Egypt (when matched to Greece) fatty acyl glycosides of mono-and disaccharides and carboxylic acids were found to increase. Greece compared to Italy showed increased levels of selected terpenoids, iridoid glycosides, and fatty acids. These differences, presented in detail below, are potentially attributed to the different orange tree varieties used in the case studies of the project, especially between Greece and Italy, where the majority of citrus-based honey came from, and the geographical origin. A recent study showed the effect of geographical origin on the differences concerning the chemical composition among the same types of honey [33]. In another domain, recent work on the chemical composition of juices originating from several European pear cultivars has shown that such differences are cultivar-dependent [34].

Flavonoids
The use of flavonoids as chemical markers of honey's botanical origin has been reported by several research groups [35,36]. Similar patterns have been used to identify bee pollen's respective origin [37]. In this work, attention was given to discriminate citrus honey from three distinct Mediterranean countries. This task was challenging since the majority of samples were of citrus crop origin and the margins of differentiation among them, at first reading, seemed limited.
Hesperidin (flavanone) and its aglycone form, hesperetin, are characteristic flavonoids of citrus honey. Hence, insight into the variability of concentrations among the same type of honey is of the utmost importance. The latter becomes even more significant when attempts are made to discriminate honeys based on geographical criteria. It is noteworthy that hesperetin (tR = 12.71, with 302.0790 m/z Table 1) was increased in Italian honeys, while hesperidin (tR = 10.14 with 610.1903 m/z Table 1), a glycoside of hesperetin, was comparatively increased in Greek honeys ( Figure 6). Both molecules can be further investigated under a biosynthetic framework. Flavonoids' biosynthesis has been studied by many research groups, though for several of its members (i.e., hesperetin), it is still unknown. Hence, any observation of their formation and the conditions that govern them is significant. It is reported that under specific conditions the content of flavanone glycosides and their aglycones undergoes drastic changes, in order to protect plants against pathogens [38]. In the specific case, since hesperetin is a precursor of hesperidin, it seems that in the specific Greek honey samples, this conversion was favored, proceeding via the formation of the intermediate hesperetin-7-O-β-D-glucoside. Hence, the observed relative difference of hesperetin levels can be attributed to such factors. To strengthen the HRMS results, the targeted HPLC-PDA-ESI/MS chemical analysis verified the differential increase in these two molecules. More specifically, hesperetin levels in Greek citrus honey showed a mean concentration of 0.19 ± 0.04 µg/g honey, while in Italian citrus honey the respective value was 0.49 ± 0.20 µg/g honey (Table S1). On the contrary, for hesperidin, the mean concentration in Greek samples was 3.78 ± 0.51 µg/g honey, while in Italian samples was 1.12±0.16 µg/g honey. In the three Egyptian citrus honey samples, hesperetin and hesperidin mean levels were 0.14 ± 0.04 and 0.71 ± 0.17 µg/g honey, respectively (  Table S1). These concentration differences might also reflect strictly regional and varietal characteristics that should not be foreseen as a generic remark for the specific type of honey (of the country of origin). The comparison with literature reports on orange blossom honey showed that the determined levels for the majority of flavonoids are in the same order of magnitude [39,40], though certain differences are observed. More specifically, naringenin's and rutin's mean concentrations in Italian and Greek honey (for rutin in Egyptian honey as well) are higher than the ones reported in Spanish honey [39], though in the specific study, only one commercial orange blossom honey sample was analyzed. It is reported that under specific conditions the content of flavanone glycosides and their aglycones undergoes drastic changes, in order to protect plants against pathogens [38]. In the specific case, since hesperetin is a precursor of hesperidin, it seems that in the specific Greek honey samples, this conversion was favored, proceeding via the formation of the intermediate hesperetin-7-O-β-D-glucoside. Hence, the observed relative difference of hesperetin levels can be attributed to such factors. To strengthen the HRMS results, the targeted HPLC-PDA-ESI/MS chemical analysis verified the differential increase in these two molecules. More specifically, hesperetin levels in Greek citrus honey showed a mean concentration of 0.19 ± 0.04 µg/g honey, while in Italian citrus honey the respective value was 0.49 ± 0.20 µg/g honey (Table S1). On the contrary, for hesperidin, the mean concentration in Greek samples was 3.78 ± 0.51 µg/g honey, while in Italian samples was 1.12 ± 0.16 µg/g honey. In the three Egyptian citrus honey samples, hesperetin and hesperidin mean levels were 0.14 ± 0.04 and 0.71 ± 0.17 µg/g honey, respectively (  Table S1). These concentration differences might also reflect strictly regional and varietal characteristics that should not be foreseen as a generic remark for the specific type of honey (of the country of origin). The comparison with literature reports on orange blossom honey showed that the determined levels for the majority of flavonoids are in the same order of magnitude [39,40], though certain differences are observed. More specifically, naringenin's and rutin's mean concentrations in Italian and Greek honey (for rutin in Egyptian honey as well) are higher than the ones reported in Spanish honey [39], though in the specific study, only one commercial orange blossom honey sample was analyzed.
Dihydrokaempferol (tR = 9.73, with 288.0634 m/z Table 1) increased in Italian honey is an important finding considering that its only report in honey is on stingless bee honey in the form of aromadendrin (its trans isomer) [41]. Hence, it can be further studied as a putative biomarker of Italian citrus honey. Plant sources for this compound mentioned in the literature belong to the families of Cactaceae and Aizoaceae [42]. The Cactacae family is also reported as flora of the case studies in Italy, flowering from April to July (see Tables S4 and S5A,B). Quercetin 3-O-sophoroside (tR = 8.15, with 626.1491 m/z Table 1) increased in Greek samples has been found at significant amounts in rapeseed honey [40] and trace levels in Diplotaxis honey [43]. Brassicaceae family of plants, from which rapeseed honey is produced is also reported in the Greek case study fields (Brassica/Sinapis, and Draba/Eruca, see also the melissopalynological analysis section and the respective Table). The same compound was also identified in Lithuanian bee pollen [44]. However, until the present work, no connection was made to citrus honey. Therefore, its use as a chemical marker for Greek citrus honey can be further studied.
Flavonoids not previously reported in citrus honey are also an intriguing task. 3,7-Di-O-methyl quercetin (tR = 13.92, with 330.0742 m/z Table 1) is a methylated flavanol and polyphenolic metabolite of Rhamnus disperma [45] (Rhamnus genus was identified via melissopalynological analysis in Italian samples, less in Greek and Egyptian samples) also reported as a constituent of Loranthus globosus, exhibiting anti-Alzheimer activity [46]. A rare flavonoid, glycoside isorhamnetin 3-rutinoside 4"-rhamnoside (tR = 10.55, with 770.2267 m/z Table 1), was putatively annotated and differentially increased in Italian citrus honey. This chemical was reported as a component of the flowers of Cucurbita pepo [47]. In general, Cucurbitaceae (pollen grains identified both in Italian and Greek honey) are plants visited by bees [48], and their fruit set is dependent on pollination [49]; hence, such chemicals are expected in apiculture commodities. Hispidulin (tR = 13.54, with 300.0635 m/z Table 1, see structure in Figure 6) is a bioactive [50] methylated-oxygenated flavone-documented component of tropical propolis [51]. To our knowledge, it is its first report on honey and apiculture commodities of non-tropical origin. Concerning other common flavonoids among the countries, their concentrations are presented in Table S1.

Fatty and Organic Acids
Interestingly, in Greek honey samples, the saturated fatty acid, 10-hydroxy-2-decenoic acid (10-HDA) (tR = 13.04, with 186.1248 m/z Table 1), a known bioactive molecule [52] was distinctively increased. More specifically, 10-HDA is the driving force behind the pronounced biological activity of royal jelly that, to an extent, justifies its increased demand from consumers. The biosynthesis of 10-HDA has gathered large attention from the scientific community. Brown and colleagues suggested sucrose as the precursor of 10-HDA after its administration in worker bees and subsequent monitoring of such fatty acids [53]. Plettner and coworkers moved forward, suggesting a de novo fatty acids synthesis with stearic acid being the precursor [54].
The differentiation in the levels of 10-HDA is an intriguing issue not only for this molecule but also for all compounds annotated in this study. Another reading for some of the constituents that are found in other apiculture matrices (such as propolis) than honey is that these bioactive compounds are transferred to honey at the onset of nectar's deposition in the honeycomb [43]. Another hydroxydecenoic acid counterpart, putatively annotated in citrus honey of the present work, is 14-hydroxy-12-tetradecenoic acid (tR = 30.19, with 242.1882 m/z Table 1). Considering the already acknowledged bioactivity of 10-HDA, the disclosure of this molecule can be an added value for Mediterranean citrus honey. (3S, 7R)-iso-jasmonic acid (tR = 13.11, with 210.1251 m/z Table 1) was differentially increased in Greek citrus honey. The latter is a signal molecule involved in the octadecanoic pathway, produced in plants, especially after insects' attacks [55]. In this context, jasmonic acid is a phytohormone reported in monofloral raw honeys [22,56], which along with other jasmonates, demonstrate anticancer activities [57]. Phytohormones are found in plants' nectar [58] and, after foraging, can be transferred inside the beehives. Ferulic acid (tR = 9.30, with 194.0572 m/z, Table 1) is a known component of honey and apiculture matrices, verified in this study. Its quantitative determination showed higher levels in Italian honey (mean levels at 1.22 ± 0.28 µg/g) compared to Greek and Egyptian ones (0.96 ± 0.21 µg/g and 0.65 ± 0.13 µg/g, respectively).
Corchorifatty acid F (tR = 14.51, with 328.2253 m/z Table 1) is a challenging molecule when viewed from the perspective of the components of honey. More specifically, this acid was reported as a constituent of the leaves of Corchorus olitorius L. (Tiliaceae) [59], and of Coryphantha macromeris (Cactaceae) [60] native in the United States of America and Mexico. Corchorifatty acid F connects with the plant Euphorbia hirta L. The latter is native to tropical and subtropical regions of Central South America, Asia, and Africa. Nevertheless, in the Mediterranean region, plants of both the Euphorbiaceae and Cactaceae families exist and recognized in the specific study (see Tables S4 and S5A,B); therefore, its putative annotation can be attributed to such plants. Last but not least, another corchorifatty acid (B) was identified in the Solanaceae family and Solanum americanum Mill., in particular [61], a family of plants abundant in Europe as well.

Terpenes
All Greek and some of the Italian honeys displayed several terpenoids in their chemical arsenal. More specifically, the qualitative analysis and semi-quantitative approach based on the relative abundances obtained (peak areas) for Greek citrus honey showed a constant representation of certain terpenoids. The latter might be attributed to the "homogeneity" and the non-scattered profile of the Greek case study fields, which are all located within an overall perimeter of 15 km. Among the compounds annotated are hallactone B (tR = 9.69, with 440.1160 m/z Table 1) and provincialin (tR = 12.68, with 518.2153 m/z Table 1, see structure in Figure 7). These terpenoids belong to the class of limonoids, which are highly oxygenated triterpenic molecules distributed in several plants' families, exemplified by the citrus crops and plants of the Cucurbitaceae families. To our knowledge for provincialin and hallactone B, it is their first report in honey. Elaborating on terpenes, it is worth noting that different terpenes were found to be elevated in Italian honey. The latter can be a point to distinct orange honey among regions and countries. However, additional sampling and data are needed to corroborate such conclusions since it is impossible to "avoid" chemical substances derived from other plants. A rare terpenic molecule among them is secologanate (tR = 6.59, with 374.1214 m/z Table 1). This compound reference is focused on the Dendrobium officinale, which is a special orchid species native to Asia [62]. Nevertheless, since Dendrobium species are cultivated in Europe, used as ornamental plants, and visited by bees, it is possible to identify them as potential sources of secologanate. Cincassiol B (tR = 12.33, with 422.1919 m/z Table 1), a diterpenic glycoside, has been reported as a constituent in Descurainia sophia seeds extract [63]. It is native to Eurasia (including Italy) but, to our knowledge, has never been reported as a honey component.
Largely different levels of the same compounds or substantially different compositions of terpenoids among crops or plants of the same type have been reported by several groups. Indicatively, terpenoids showed dissimilar contents among different types of tangerines [64]. In another study, the terpenoid variations were also correlated to genotypic differences [65]. Consequently, the results of the differential expression of terpenoids in the present study are not unexpected.

Other Chemicals
Lumichrome (tR = 10.86, with 242.0802 m/z Table 1) is the enzymatic cleavage product of riboflavin. In this context, its annotation in the present work comes to add another evidence to its reference to monofloral honeys [66]. The same group used this molecule as a biomarker of thistle honey [67]. Though formononetin (tR = 14.15, with 268.0737 m/z Table 1) has been ascribed to acacia, vitex, and linden honey from China [62], this report is its first appearance in citrus honey. Formononetin has also been reported as a constituent of Brazilian red propolis [68], a natural bee product known for its antibacterial, antioxidant, and cytotoxic activities. A precursor of sphingolipids known as parts of biological membranes, sphingosine, was also identified in the present study. The major phospholipid of the honeybee is sphingomyelin, reported since the early 1970s [69]. Greek honey differentially increased (4E, 6E, d14:2) sphingosine (or tetradecasphinga-4E,6E-dienine), (tR = 30.19, with 241.2041 m/z Table 1)), an amino alcohol first reported in citrus honey in this study.
4-Hydroxycinnamyl aldehyde (tR = 9.06, with 148.0514 m/z Table 1) was differentially increased in Italian citrus honey. The latter is a cinnamaldehyde derivative isolated from the rhizomes Alpinia galanga (Linn) [70] and in Salix species [71]. Salix alba L. was reported in the flora of the Italian fields; hence, it can be the potential source of this molecule, though not. It is noteworthy that it is a bioactive molecule inducing the apoptosis of human leukemic cells [72]. Scopoletin (tR = 8.09, with 192.0424 m/z Table 1), which is reported as a constituent of Greek cotton honey [73], has not been evidenced in citrus honey until the specific report. In the same context, cotton cultivation is not registered in the Greek case study areas to imply the uptake of this component from such an antagonizing crop. However, it is found in Asteraceae, Convolvulaceae, Rubiaceae, Solanaceae, and Moraceae plants (see [74] and references therein), with Asteraceae, Convolvulaceae, and Moraceae, verified in this work. Other components, such as sebacic acid (tR = 11.31, with 202.1206 m/z Table 1), are fermentation products of the bee gut microbiota [75]. Acetophenone (tR = 10.81, with 120.0576 m/z Table 1) belongs to the volatile markers of several types of honey, such as chestnut [76], thyme [77], raspberry, heather, and rape honey [78]. Consequently, its relatively high abundance in Italian citrus honey is worth pointing out, though further use as a chemical marker seems limited. Abscisic acid (tR = 11.07, with 264.1362 m/z Table 1) was relatively increased in Italian citrus honey, exhibiting a mean concentration of 10.01 ± 0.45 µg/g honey, compared to 8.78 ± 0.39 µg/g and 6.86 ± 0.49 µg/g in Greek and Egyptian honey correspondingly. The latter is a known plant hormone and chemical marker of heather honey [22,79].
4-Hydroxyquinoline (tR = 4.78, with 145.0527 m/z Table 1) increased in Greek honey samples expands the heterocyclic portfolio of citrus honey. Quinolines are bioactive components (especially their fused alkaloid structures) rarely reported in honey [80]. However, it is imperative to mention that the specific compound was used as a chemical marker of jujube honey [81]. In addition, it is reported as a constituent of Rutaceae plants, abundant in Greek case study areas.
Molecules 2023, 28, x FOR PEER REVIEW though not. It is noteworthy that it is a bioactive molecule inducing the a man leukemic cells [72]. Scopoletin (tR = 8.09, with 192.0424 m/z Table 1), wh as a constituent of Greek cotton honey [73], has not been evidenced in citr the specific report. In the same context, cotton cultivation is not registere case study areas to imply the uptake of this component from such an anta However, it is found in Asteraceae, Convolvulaceae, Rubiaceae, Solanace ceae plants (see [74] and references therein), with Asteraceae, Convolvulac ceae, verified in this work. Other components, such as sebacic acid (tR 202.1206 m/z Table 1), are fermentation products of the bee gut microbiota none (tR = 10.81, with 120.0576 m/z Table 1) belongs to the volatile markers o of honey, such as chestnut [76], thyme [77], raspberry, heather, and rape h sequently, its relatively high abundance in Italian citrus honey is worth though further use as a chemical marker seems limited. Abscisic acid (tR 264.1362 m/z Table 1) was relatively increased in Italian citrus honey, exh concentration of 10.01 ± 0.45 µg/g honey, compared to 8.78 ± 0.39 µg/g and in Greek and Egyptian honey correspondingly. The latter is a known plan chemical marker of heather honey [22,79].
4-Hydroxyquinoline (tR = 4.78, with 145.0527 m/z Table 1) increased i samples expands the heterocyclic portfolio of citrus honey. Quinolines are ponents (especially their fused alkaloid structures) rarely reported in ho ever, it is imperative to mention that the specific compound was used as a ch of jujube honey [81]. In addition, it is reported as a constituent of Rutacea dant in Greek case study areas. In the same framework, iridoid components identified increased in th honey are reported for the first time in the present study. Iridoids are glyco various plants, and they bind to glucose. They have the general form of cy and a molecular structure related to iridodial. Nepetaside (tR = 11.87, wit Table 1) has been reported as a constituent of Gentiana [82] and Nepeta sp these plant species are reported in Greece [84]. Recently an RP-HPLC-MS nepetaside as a secondary metabolite of Heliotropium crispum Desf. [85], w spread genus of the Boraginaceae plant family in the Mediterranean region In the same framework, iridoid components identified increased in the Greek citrus honey are reported for the first time in the present study. Iridoids are glycosides found in various plants, and they bind to glucose. They have the general form of cyclopentopyran and a molecular structure related to iridodial. Nepetaside (tR = 11.87, with 346.1629 m/z Table 1) has been reported as a constituent of Gentiana [82] and Nepeta species [83]. Both these plant species are reported in Greece [84]. Recently an RP-HPLC-MS study unveiled nepetaside as a secondary metabolite of Heliotropium crispum Desf. [85], which is a widespread genus of the Boraginaceae plant family in the Mediterranean region. Therefore, the occurrence of this rare iridoid can be attributed to plants of this family that exhibited high frequency in Greek honey samples. For patrinoside (tR = 10.37, with 508.2155 m/z Table 1), another iridoid glycoside, a known component of Valeriana species [86], the picture is clearer. Hence, its annotation in honey samples from Greece can be substantiated based on the occurrence of plant species such as Sambucus ebulus L. (Adoxaceae, also found in forestall areas along with fir tree) and Centranthus sp. (abundant in the Argolis region, based also on optical observation) that contain this component and its aglycone forms [87][88][89]. Other plants that contain this molecule [90,91] are not reported in the Greek flora. Like terpenes, a different iridoid was distinctively increased in Italian honey than the Greek ones. The detection of iridoids, an underexplored chemical family in honey, is of utmost importance, considering their biological properties and the added value that they induce in the natural products containing them [92]. Other compounds associated with citrus honey, flowers, and nectar or citrus juices [93], such as caffeine [94], exemplified by synephrine (a biogenic amine) also used as a selective marker of citrus honey [95], were not detected in the presented work via LC-HRMS, despite their inclusion in the HPLC-PDA-ESI/MS targeted analytical method. Consequently, the presented work substantiates the complexity of secondary metabolite generation among the same matrices and their use as exclusive chemical markers.
The majority of plant species reported in this work as potential sources of chemicals annotated, are visited, and pollinated by bees for nectar and pollen collection, therefore the occurrence of these plants' components in bee products is logical. To substantiate this statement, Gao and colleagues investigated the bioactive arsenal of citrus nectar and honey [96]. Compounds detected in nectar, such as hesperetin, rutin and gallic acid (see Tables 1 and S1) were also identified in the present study corroborating a high likelihood of their transfer from citrus nectar to honey. Last but not least, flower-plant traits are associated with the motifs of bees' visits and plant-pollinator interactions [97,98]. Specifically, the variations and frequency of bees' visitation (and population consequently) in the diversity of flowers that come across during foraging, can affect the degree of chemicals transferred to honey. Hence, it is an additional factor shaping honey's chemical profile, apart from nectar and pollen composition.

GC-MS Findings
GC-MS analytical results (no full data shown) were in line with the literature findings on citrus honey volatile components. More specifically, methyl anthranilate, lilac aldehydes (all isomers I-IV), and 2,6-dimethyl-2,7-octadiene-1,6-diol, all characteristic markers of citrus honeys [13,14] were identified in this study as well. Similarly, a plethora of saturated fatty acids and common constituents of citrus honey was identified, such as dodecanoic, decanoic, and undecanoic acid, and other aldehydes such as hexadecanal and heptadecanal. Linalool was also putatively identified in some honey samples. Nevertheless, some volatile chemicals firstly reported in citrus honey are worth mentioning. More specifically, and to our knowledge, 5-dodecyldihydro-(3H)-furanone and 5-tetradecyldihydro-2(3H)-furanone (furan derivatives) were identified for the first time in Greek citrus honeys. Dihydrofuranones are reported in several types of honey [76], but not the aforementioned.

Melissopalynological Analysis
The relative level of abundance and relative frequency of citrus pollen in unifloral citrus honeys classifies it as under-represented pollen [12,99]. This means that a relatively small percentage of citrus pollen grains is expected in the melissopalynological analysis, with the remainder of pollen representing various plants of the local flora. The legal limits for the relative frequency of citrus pollen in the monofloral citrus honey range between 3% and 20%, regarding the national legislation of several European countries [100].

Melissopalynological Analysis of Citrus Honey from Sicily, Italy
In citrus honey from Sicily (n = 20), 56 plant families were identified (33 of nectariferous and 23 of nectarless plants, Table S5A). The families of which the sum of pollen accounted for 95-100% of the whole pollen spectra of nectariferous plants are shown in Table 2. The Rutaceae family (Citrus) was found, as expected, in low relative frequencies (m.v. 6%, range: <1-31%, Table 2). The dominant family was the Fabaceae, present in all samples examined (n = 19, m.v. 49%, range 15-92%, Table 2). Ten genera of this family were identified. Among them, Lotus and Melilotus were predominant (>45% in some samples), while Trifolium repens type and Vicia were secondary (Tables 2 and S5A). Other important families were the Asteraceae, Boraginaceae, Brassicaceae, and Rosaceae. Notable was the Boraginaceae family, observed in 15 out of 19 honey samples. Five plant genera of this family were identified, and among them, Echium and Cynoglossum/Cerinthe were classified as secondary pollen (16-45%) in some samples (Tables 2 and S5A). Other genera with secondary pollen (at least in one sample) were Sinapis/Brassica and Draba/Eruca (Brassicaceae), Cirsium and Centaurea solstitialis type (Asteraceae), Ferula-Scandix (Apiaceae). All other plant species found in citrus honey samples were in the range of 3-16% (important minor pollen) or <3% (minor pollen) (Table S5A,B).
In the quantitative melissopalynological analysis, 75% of the samples provided results compatible with the class of unifloral honeys with under-represented pollen (Class I: N ≤ 20 × 10 3 ), while the rest of the samples were measured in Class II (N = 21 × 10 3 -100 × 10 3 ), where N is the number of plant elements in 10 g honey (Table S5B). In the specific samples, N is in fact the number of pollen grains, as no significant amount of honeydew elements was detected [99]. Finally, the honey samples from the Italian market (n = 12) did not differ significantly from the Italian honeycomb honeys (n = 8), except for one that was excluded from statistical analysis ( Figure 8). Nevertheless, a great variety of plant species was observed in citrus honeys from the fields (Tables 2 and S5A,B).

Melissopalynological Analysis of Citrus Honey from Argos, Greece
In citrus honey from Argos, Greece (n = 12), the melissopalynological analysis revealed 51 plant families (28 of nectariferous and 23 of nectarless plants, Table S5A). The families of which the sum of pollen accounted for 96-100% of the whole pollen spectra of nectariferous plants, are shown in Table 3.
Most of the families are common to those encountered in citrus honeys from Italy. However, notable differences exist, such as the relative abundance of each family and the number of plant species per family. The pollen of genus Citrus is considered underrepresented in honey [12,99]. Nevertheless, in citrus honeys from Argos, it was classified as predominant pollen (m.v. 53%, range 17-79%, Table 3). This is mainly due to the absence of other nectariferous plants of the same flowering period as the orange trees, but the variety of orange trees in the field can also interplay. Another fact that affects the relative frequencies is the large percentage of pollen of nectarless plants (m.v. 84%, range 68-91%, Table S5A,B), which is excluded during the calculation of the frequencies of nectariferous plants. Apart from the dominant Rutaceae family (Citrus), other important families were the Brassicaceae, Asteraceae, and Boraginaceae families ( Table 3).
The genus Sinapis/Brassica (Brassicaceae) was classified as secondary in one-third of the samples, and the genus Echium (the only representative of the Boraginaceae family in Greek citrus honey samples) was secondary or even dominant in some samples. Anthemis and Taraxacum pollen types were the most abundant in the Asteraceae family. Other plant species with a high relative frequency (dominant or secondary, at least in one sample) were the Pyrus-Prunus type (Rosaceae) and Tordyllium type (Apiaceae) ( Tables 3 and S5A).
Among the 23 families of nectarless plants encountered in Greek citrus honeys, the most important were the following (in brackets the corresponding genera): Fagaceae (Quercus coccifera type), Oleaceae (Olea), Anacardiaceae (Pistacia), and Papaveraceae. The genera Quercus and Olea were dominant (>45% of all pollen types) in several samples. The average percentage of nectarless plants was high (m.v. 84%, range 68-91%, Table S5A,B). The whole pollen spectrum, observed in all honey samples (n = 12), was consistent with the local flora [102,103]. In the quantitative melissopalynological analysis, only one sample was classified as Class I (≤20 × 10 3 ), nine samples as Class II (21 × 10 3 -100 × 10 3 ), and one sample was classified as Class III (101 × 10 3 -500 × 10 3 ) [99]. The under-represented nature of citrus pollen could be confirmed only if the pollen of nectariferous plants was taken into account (see PG-nectariferous/10 g, m.v. 6145, range 2480-9118, Table S5B). This is due to a large number of pollen grains of nectarless plants that do not actually contribute to the formation of honey, but they come from pollen already stored in the honeycomb by the bees. Finally, the relative frequencies of pollen, measured in the seven honey samples and the five honeycomb samples, were quite different, indicating that the presence of honeycomb can affect to a great extent the results of melissopalynological analysis (Tables S5A,B and 3).

Melissopalynological Analysis of Citrus Honey from Al Shaeir Island, Egypt
In honeys from Egypt (n = 8), the melissopalynological analysis revealed the existence of 34 families (19 of nectariferous and 15 of nectarless plants, Table S5A). Although five honeys were declared as unifloral (2 citrus honeys, 2 clover honeys, and 1 basil honey) and three honeys as polyfloral; mainly from aromatic medicinal plants and citrus trees, the most important families were common to all samples (Table 4).
Despite the under-represented nature of Citrus pollen, its relative abundance in the two citrus honeys was high (79% and 29%, respectively, Table 4). In the two clover honeys, the relative frequency of Trifolium sp. pollen was 83% and 25% and finally the relative frequency of Ocinum pollen in the basil honey was <1% (Table 4). This could be explained by the fact that basil flowers, although of great apicultural value, display a poor pollen collection by bees [104].
The most important families, of which the sum of pollen accounted for 95-100% of the whole pollen spectra of nectariferous plants, are shown in Table 4. Apart from Fabaceae and Rutaceae family, expected to be observed in high relative frequencies, as Trifolium sp. and Citrus trees are cultivated in the specific fields, notable was the Myrtaceae family. In half of the samples, Eucalyptus pollen was dominant (>45%), and in two samples was secondary pollen (16-45%). From the Brassicaceae family, Lepidium/Draba type was dominant pollen in one sample and secondary in another. The Malvaceae family, although in low relative frequencies (<2%), was observed in 6 out of 8 samples, represented by more than three plant genera (Bombax, Abutilon, Sida, etc.). The only representative of Lamiaceae family was genus Ocinum, encountered in several samples, though in low frequencies (<1%). Finally, the Asparagaceae family was observed in one sample as dominant pollen (Tables S5A,B and 4).  Among the 17 families of nectarless plants, the most important were the following: Poaceae (Graminae), Arecaceae (Palmae), Solanaceae, Cyperaceae, Asteraceae, Chenopodiaceae/Amaranthaceae, Cistaceae, and Oleaceae. The Arecaceae (Palmae) family was observed as secondary or important minor pollen in several samples and the Asteraceae family (Artemisia and Xanthium genus) was dominant pollen in the basil honey. The average percentage of nectarless plants was relatively low (m.v. 23%, range 1-64% of all pollen types, Table S5A,B). The whole pollen spectrum, observed in all honey samples (n = 8), was consistent with the local flora and the plants cultivated in the specific fields [105].
Notable differences from the Italian and Greek honeys were the strong presence of the Myrtaceae family, the various genera of the Malvaceae family and the predominance of nectarless plant families different from those of the Italian and Greek samples, especially the Arecaceae (Palmae) family (Table S5A).
Quantitative melissopalynological analysis provided contradictory results, as a large variation was observed: one sample was classified in Class I (N ≤ 20 × 10 3 ), two samples in Class II (N = 21 × 10 3 -100 × 10 3 ), three samples in Class III (N = 101 × 10 3 -500 × 10 3 ) and two samples in Class V (N > 10 6 ). A possible explanation for this, is that honey was extracted from the honeycomb by pressing or squeezing, thus making possible the enrichment of honey with pollen stored in the honeycomb by the bees, possibly in cells blocked by a quantity of honey.
Overall, citrus honeys from the three different countries (Italy, Greece, and Egypt), as well as honeys of different botanical origins from Egypt, showed a different melissopalynological profile, which made it possible to distinguish the samples using principal component analysis (Figures 8 and 9). With the exception of one sample from Italy (CH-IT-20) and the "basil honey" from Egypt, all other samples were well grouped.

Environmental Factors
The biosynthesis and accumulation of plant secondary metabolites (PSMs) a erned by inherent and extrinsic factors. The latter encompass biotic and abiotic e mental drivers such as the response to parasites, insects, light irradiation, tempe water content on soil, etc. [106].
Research has shown that the selectivity of biosynthetic reaction pathways in p influenced by seasonal variations reporting differences in flavonoid levels and [107]. Despite not the scope of the present research, environmental parameters we templated in this work in an effort to explain the obtained results.
Annual mean temperatures, humidity, rainfall, sunlight, etc., were acquired fr nearest meteorological station or from the National Meteorological Data Service Ce Greece, Italy, and Egypt. The latter provides the most precise data for the respec gions. Α long period of light irradiation was reported to increase the concentrat PSMs, such as flavonoids and phenolic acids. Light intensity is also a parameter disregard.

Environmental Factors
The biosynthesis and accumulation of plant secondary metabolites (PSMs) are governed by inherent and extrinsic factors. The latter encompass biotic and abiotic environmental drivers such as the response to parasites, insects, light irradiation, temperature, water content on soil, etc. [106].
Research has shown that the selectivity of biosynthetic reaction pathways in plants is influenced by seasonal variations reporting differences in flavonoid levels and content [107].
Despite not the scope of the present research, environmental parameters were contemplated in this work in an effort to explain the obtained results.
Annual mean temperatures, humidity, rainfall, sunlight, etc., were acquired from the nearest meteorological station or from the National Meteorological Data Service Centre of Greece, Italy, and Egypt. The latter provides the most precise data for the respective regions. A long period of light irradiation was reported to increase the concentrations of PSMs, such as flavonoids and phenolic acids. Light intensity is also a parameter not to disregard.
Nevertheless, the three Mediterranean regions share common weather features, adding complexity to the challenging task of attributing specific chemicals' relative increase or decrease to environmental factors. For the Argolida region in Greece, an average temperature of 15.4 • C and relative humidity of 67.1% was recorded for March-May 2020 (total annual sunshine reaching 2561 h per year was also documented). The solar radiation observed in the area (1700-1800 kWh/m 2 /y) is higher in Greece. In Sicily (Csa, hot summer Mediterranean climate), where 8 experimental plots were established, the annual sunshine reaches 2700 h that is 1.05 times higher than the respective in Argolida Greece. The respective regions in Egypt are similarly characterized by a hot desert climate (BWh) which is extremely dry with virtually no rainfall. Annual sunshine in Egypt is higher than in Greece and Italy, with 3300 h in the upper part (close to the Mediterranean). More specifically, in Egyptian pilot areas the period March to May 2020 an average temperature of 21.7 • C, and a relative humidity of 49.3% were recorded. Based on the literature long photoperiod increases the biosynthesis and concentrations of flavonoids and phenolic acids, compared to short durations of light irradiation [106]. Flavonoids and phenolic acids were less expressed than in Italy and Greece. However, for Egypt, the latter was not verified, a fact than can also be attributed to the lower number of Egyptian citrus honey integrated into the specific study, and the less diverse flora of the Egyptian pilot areas. On the other hand, the metabolomics approach to Greek and Italian honey might postulate that apart from the citrus tree varieties, the higher average sunshine in Sicily might have favored the formation of more flavonoids and phenolic acids than in Greece, considering that both areas were rich in AMPs. On the contrary, terpenic molecules are constantly differentially increased in Greek honey. The latter can possibly be attributed to the more shaded conditions and lesser sunshine that preponderate in the region compared to Italian region, which is reported to favor the formation of some terpenes [108]. Nevertheless, such a conclusion needs careful consideration since it cannot be regarded as a generic rule for all terpenes since studies showed that not all terpenes biosynthesis are affected the same way under comparable environmental conditions [109]. In addition, it reflects the specific honey samples collected during 2020, whose chemical portfolio through the biosynthesis of all classes of plant secondary metabolites (terpenes, phenolic and nitrogen-containing compounds) is shaped predominantly by the citrus crop, the flora of the region and the climatic conditions.
Since the case study areas are not homologous, though most of them refer to the predominating citrus crop and derived citrus honey, it is risky to further elaborate on how other parameters (solar radiation in connection to altitude) affected the generation of PSMs. Nevertheless, this work adds another piece of evidence on the multifactorial issue of PSMs biosynthesis and the variations in the diversity of compounds and their concentrations.
A deeper analysis of individual flavonoids, phenolic acids, terpenoids, and other chemicals' concentrations and variations will continue to further reinforce the presented findings. To build upon this work, additional honey samples from the two succeeding years (2021 and 2022, beehives were also placed in the same case study areas) of the PLANT-B project will be incorporated, and the impact of installed AMPs might be more evidenttraceable in honey's composition. The latter will be viewed under an overall metabolomics scheme along with prospective and solid quantification of not yet identified or reported bioactive compounds.

Honey Samples
A total of 28 samples were collected from all countries: Greece (7 honey + 5 replicate honeycomb honey), Italy (8 honeycomb honey), and Egypt (8 honey). To reinforce the metabolomics approach, 12 citrus honey samples from the Italian market were purchased and included in the study. The majority of samples were citrus honey, yet some honey samples were produced exclusively from AMPs orchards. More specifically, in Egypt, two fields featured citrus in addition to clover in one of them, and the second field had clover in addition to the combination of basil, borage, coriander, anise, and caraway. In Greece and Italy, all samples were derived from citrus orchards, some of them containing installed AMPs.
In Greece, 7 experimental orange orchards, cv Navelina and Merlin, were established in the region of Argolida (see Figure S1). In each field, five 10-frame Langstroth hives beehives were placed. The beehives were placed on pallets. The bee colony hives were of the same honeybee population, with queens of the same age, and had been handled in the same way during the winter and spring. The installation of the beehives in the experimental fields occurred on 23 April 2020, while the orange trees were at the beginning of flowering (5-10%) in most orchards and in 30-50% of full flowering in orchards of early blooming. The beehives were inspected every week, and the necessary beekeeping manipulations were carried out (addition of frames, swarming control, expanding space for honey and egg laying, etc.) and were recorded. The honey harvest was made on 14 May, 3 weeks later and while the orange trees were in 50-100% of full flowering.
The frames with the ripe honey placed on empty floors were transferred on the same day to the Laboratory of Apiculture of the Institute of Mediterranean Forest Ecosystems and Forest Products Technology (Greece), and the next day the extraction of honey occurred. Two kinds of samples were collected, honey samples from the honey extract and honeycomb samples with honey. The pieces of honeycomb with sealed honey were placed in glass jars. Honey from the honeycombs was obtained at Benaki Phytopathological Institute (Greece) by squeezing, and the honey was strained by gravity into a clean container.
In the frames that were placed in the honey extractor, wax cappings were removed with a cold scratcher, and the honey was also strained by gravity from a stainless strainer to a honey tank. The honey was then allowed to rest for at least one week, during which time all air bubbles present in the honey floated to the top. Then pure honey was drained through a latch at the bottom of the tank to fill glass jars. All samples were placed in a freezer at −18 • C until the day of their analysis.
In Italy, the Sicilian experimental orange orchards ( Figure S1) contained the Navel orange variety, with the exception of the experimental sites in the municipality of Chiaramonte Gulfi (IPM3, O1, and O2) that contained the "Tarocco" orange varieties. Honey sampled during the first year of experimentation was derived from Italian-Dadant hives with 10 and 6 frames. The supers were removed by preparing the bee escapees 24 h before. Subsequently, the wax capping was removed, and the frames were centrifuged. The extracted honey was filtered and decanted. Occasionally, some hives had not filled the supers, and therefore, the honey was recovered directly from a portion of the honeycomb. Egyptian experimental orange orchards ( Figure S1) also contained Navel oranges. The complete information on samples, location, and installed AMPs are presented in Figure S1 and Table S1 (Supplementary Material).

Sample Preparation
Honey (3 g) was homogenized with the aid of a glass rod in 7 mL of acidified ultrapure H 2 O (pH = 2.2). Then, the mixture was loaded on a pre-activated DSC-C18 SPE cartridge (activation was performed by succeeding elution of MeOH (3 mL) and H 2 O (3 mL)). Consequently, the SPE cartridge was washed with acidified H 2 O (2 mL) and 5 mL of ultrapure H 2 O. Compounds of interest were eluted with 3.5 mL of a MeOH:ACN (2:1, v/v) solution. Then, the eluate was evaporated to dryness using a rotary evaporator (bath temperature not surpassing 30 • C) and reconstituted in MeOH:H 2 O (7:3, v/v). After filtration with nylon filters (0.22 µm), the extract was injected into LC-HRMS. A pooled quality control (QC) sample was prepared by transferring 10 µL of each sample to an LC vial to assess system stability throughout the batch analysis.

Ultra-High-Performance Liquid Chromatography Coupled to Orbitrap High-Resolution Mass Spectrometry Analysis of Honey Extracts
UHPLC was performed using a Dionex Ultimate 3000 UHPLC system (Thermo Scientific, Karlsruhe, Germany). A Hypersil Gold UPLC C18 (2.1 × 150 mm, 1.9 µm) reversed phased column (Thermo Scientific, Germany) was used for the separation of the analytes. The analysis was performed on a Q-Exactive Orbitrap mass spectrometer using a negative and positive heating electrospray ionization source (Thermo Scientific, Germany).
The mobile phase consisted of solvents A: aqueous 0.1% (v/v) formic acid and B: MeOH. A gradient elution methodology from 0 to 40 min has been employed as follows: T = 0 min, 20%B; T = 2 min, 20%B; T = 12 min, 70%B, T = 32 min, 95%B, T = 37 min, 95%B, T = 37.1 min, 20%B; T = 40 min, 20%B. The flow rate was 0.220 mL/min, and the injection volume was 5 µL. The column temperature was kept at 35 • C while the sample tray temperature was set at 10 • C.
The optimized conditions for analysis were set as follows: capillary temperature, 350 • C; spray voltage, 2.7 kV (for negative) and 4 kV (for positive); S-lense Rf level, 50 V; sheath gas flow, 40 arb. units; aux gas flow, 5 arb. units; aux. gas heater temperature, 50 • C. Analysis was performed using the Fourier transform mass spectrometry mode of the LTQ orbitrap (FTMS) in the full scan ion mode, applying a resolution of 70,000, while the acquisition of the mass spectra was performed in every case using the centroid mode. The mass range for full MS was set at 120-1200 m/z. The data-dependent acquisition capability has been used at 35,000 resolution, allowing for MS/MS fragmentation of the three most intense ions of every peak applying a 10 s dynamic exclusion. Stepped normalized collision energy was set at 40, 60, and 100. The Xcalibur version 4 was used for data acquisition. A pooled QC sample was analyzed three times in the beginning and three times at the end of the acquisition, and every after six samples.

High-Performance Liquid Chromatography Photo Diode Array Electrospray Mass Spectrometry (HPLC-PDA-ESI/MS)
A Shimadzu (Kyoto, Japan) LCMS-2010 EV Liquid Chromatograph Mass Spectrometer instrument was used with the LCMS solution version 3.0 software consisting of a SIL-20A prominence autosampler and an SPD-M20A diode array detector. The latter were coupled in series with a mass selective detector equipped with an atmospheric pressure ionization. The LC separation was achieved on a Zorbax Eclipse Plus, 3.5 µm, 150 × 2.6 mm i.d. chromatographic column (Agilent Technologies, Santa Clara, CA, USA). The mobile phase consisted of two channels, channel: 0.1% formic acid in water (A) and ACN (B). The flow rate and mobile phase gradient were identical to previously published work [110]. Similarly, the validation of the analytical method and quantitation of analytes was grounded on the previous work of our group [110], adapting the same sample preparation followed for the UHPLC-HRMS analysis.

GC-MS Analysis SPME Holder and Fibers
The SPME holder and coated fibers (85 µm PolyAcrylate (PA), 100 µm Polydimethylsiloxane (PDMS), and 65 µm CarboWax/divinylbenzene (CW/DVB)) were supplied by Supelco (Bellefonte, PA, USA). For honey extraction, PA fiber was selected. More specifically, 1 g of honey was mixed with 5 mL 10% NaCl, shaken in a vortex mixer for 1 min, placed into a 9 mL headspace vial (Alltech, Alltech Ass. Deerfield, IL, USA), and positioned in an aluminum block heater. After 30 min of the preheating period at 60 • C with simultaneous stirring, the SPME needle penetrated the vial septum, and the fiber was exposed in the headspace of the solution. Sampling was performed for an additional 20 min at 100 • C, and finally, the needle was removed from the vial and injected into the heated injection port of the gas chromatograph. A Shimadzu Nexis GC 2030 gas chromatograph (Shimadzu Corporation, Kyoto, Japan) equipped with an AOC-6000 autosampler and a Shimadzu GCMS-TQ8040 NX triple quadrupole. Data acquisition and processing were performed by LabSolutions GCMS solution software, version 4.52. Desorption was allowed for 5 min. The injector port temperature was kept at 250 • C. A MEGA 5-HT (MEGA S.r.l., Legnano, Italy) column (30 m length × 0.25 mm i.d. × 0.25 µm film thickness) was used. The instrument worked at a constant flow of 1.55 mL/min, using the full scan mode (range of 50-500 amu) and helium (99.999% purity) as the carrier gas. The GC oven program was the following: Initial temperature at 40 • C (stayed for 3 min), ramped linearly to 260 • C, at a rate of 5 • C/min, where it stayed for an additional 8 min (overall runtime 55 min).

Data Post-Processing and Chemometrics Analysis
Raw data were exported to Compound Discoverer 2.1 (Thermo Fisher Scientific Inc., San Jose, CA, USA) for peak detection, deconvolution, deisotoping, alignment, gap filling, and composition prediction procedures. Peaks in which the quality control (QC) sample coverage was less than 50% and the relative standard deviation of the areas under the peaks was more than 30% were excluded. The normalization was QC-based, applying the median absolute deviation (MAD) normalization type. Adjusted p-values were calculated using the Benjamini-Hochberg to display the statistically significant variables. Log2 fold change values were calculated to show the variation of the selected differential metabolites between groups. The generated peak lists (for positive and negative mode) from Compound Discoverer containing the accurate masses and retention time paired with corresponding intensities for all detected peaks of all the samples were exported as a .csv file, imported to Microsoft Excel 2010 and manipulated appropriately using the "concatenate", "round" and "transpose" commands. These lists were imported in R to generate heatmaps based on the geographic origin of the honey extracts, using the function heatmap.2 from the Bioconductor package, gplots.
The data were then subjected to multivariate statistical analysis (MVA), i.e., principal component analysis (PCA) and orthogonal projections to latent structures discriminant analysis (OPLS-DA) using the ropls R Bioconductor package [111] in order to determine the optimal number of components (Supplementary Table S6), confirm the validity of the model by permutation testing utilizing 200 random permutations, detect outliers and perform feature selection with variable importance in projection (VIP) scoring from OPLS-DA models. For both PCA and OPLS-DA, different scaling methods were assessed, such as mean-centering only, mean-centering with pareto scaling, and mean-centering with unit variance scaling after log10 data transformation. For both analyses, the standard scaling method was selected as it afforded better clustering between the groups. For the selection of the markers that are differentially increased among the countries, the features that exhibited a VIP scoring greater than 1.5 were further verified by adjusted p-value (≤0.05) and Log2 fold change (>1.5). These variables were characterized as the variables that contribute the most to the clustering formation observed in OPLS-DA analysis.
For the putative annotation of the compounds the mzCloud (https://www.mzcloud.org, accessed on 10 June 2021) database was used applying m/z tolerance of 5 ppm taking into consideration the isotope distribution similarity and MS/MS fragmentation pattern.
For the compounds for which there was no record, the online in silico fragmentation tool MetFrag [112] was used, applying 5 ppm search tolerance and 0.001 mass deviation to match generated fragments against MS/MS peaks. For the MetFrag workflow, the candidate structures were retrieved from the databases Kegg (http://www.genome.jp/ kegg/compound/, accessed on 10 June 2021), LipidMaps (www.lipidmaps.org, accessed on 10 June 2021), and the Human Metabolome Database (https://hmdb.ca, accessed on 10 June 2021). For the selection of the compounds, a > 0.8 final score was applied.

Melissopalynological Analysis
For the melissopalynological analysis, the methods of Louveaux and von der Ohe were applied [99,113]. In this respect, a solution of 10 g of honey in 20 mL of distilled water was centrifuged for 10 min at 1000× g. After decanting the supernatant liquid, the sediment was suspended again in 20 mL of distilled water and centrifuged for 5 min at 1000× g. The supernatant liquid was decanted, and the sediment was spread over a microscope slide on a heating plate (40 • C). The sediment is left to dry and then is colored with some drops of a solution of fuchsine (~0.05‰ w/v in ethanol/water 1/1). After drying (40 • C), a coverslip with a drop of glycerin jelly (mounting medium; Kaiser's Glycerol Gelatin TM Merck 1.09242.0100) was placed upon the sediment, and the slide was left on the heating plate for another 5 min. The microscope slide is investigated under the microscope (400×). For the determination of relative frequencies of pollen types, 500 to 1000 pollen grains were counted [114]. To detect all possible plant sources, especially those that exhibit very large pollen grains, which are not expected to be present in a substantial amount, the entire slide was scanned at 50× magnification. Quantitative melissopalynological analysis was achieved using Slide Grids, 20 mm × 20 mm. At least 10 mm 2 were counted.
The assignment of pollen phenotypes to plant genera or species was accomplished using electronic atlases of pollen (pollen databases) [104,[115][116][117][118] and related articles [105,119,120]. It was also confirmed by the project participants through their local collaborators (in Italy, Egypt, and Greece) as well as by plant databases [103,[121][122][123]. Especially for the case of Greek honeys, in-house databases were also used. In any case, the genus or species referred to in the text and tables correspond to a pollen phenotype. The multivariate analysis of the melissopalynological analysis results was conducted using Minitab ® 17.1.0. Software. A total of 136 variables were used for the 40 samples examined. A number of 39 of the new variables were able to describe 100.0% of the system.

Conclusions
The application of an LC-HRMS metabolomics workflow combined with an elaborate melissopalynological analysis managed to unveil several new potential markers of Mediterranean citrus honey potentially associated with citrus crops and the local indigenous flora. Even though most samples were of citrus origin, 40 differentially increased compounds emerged, suggesting a potential effect of both the geographical and botanical origin of the samples. Italian samples were principally characterized by relatively increased levels of flavonoids (such as hispidulin, dihydrokaempferol) and glycerophospholipids, while in Greek samples, selected terpenoids (such as provincialin), iridoid glycosides (e.g., patrinoside), and fatty acids predominated. Egyptian samples showed augmented levels of suberic acid and fatty acyl glycosides. The different varieties of citrus crops, especially between Italy and Greece, seem to shape these differences, subsidized by the local flora and the geographical origin. Concurrently, quantification of common compounds by targeted HPLC-PDA-ESI/MS was also conducted, showing comparable results with the literature. Even though not at the forefront of this work, GC-MS managed to disclose two furan derivatives, never reported, to our knowledge, in citrus honey, as distinguishing chemicals of Greek citrus honey. Consequently, the present work managed to identify molecules that can function as chemical markers of citrus honey, covering a broad area of the Mediterranean basin. Newly reported markers not only better characterize honey but can also help beekeepers to improve the position of their product in a global antagonistic market where quality demands additional scientific evidence.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/molecules28093967/s1, Figure S1: Sampling locations of the PLANT-B case studies in Greece (Argolis, Peloponnese) Egypt (Al Shaeir Island), Italy (Sicily) (Images obtained from Google Earth); Figure S2: Loading plots of Melissopalynological parameters in the space of the two first principal components Table S1A: Egyptian case study honey samples; Table S1B: Italian case study honey samples; Table S1C: Greek case study honey samples and exact location; Table S2: Flavonoids, phenolic and organic acids quantified by HPLC-PDA-ESI/MS in orange blossom honey from the three countries; Table S3: The dominating flora during the citrus flowering season of experimental fields at Al-Qanater Al-Khairiya station-Alshaeir Island-Egypt [124,125]; Table S4: The flora of experimental fields in Sicily, Italy; Table S5A: Qualitative melissopalynological analysis for Greece (EL), Italy (IT), and Egypt (EG); Table S5B: Quantitative melissopalynological analysis for Greece (EL), Italy (IT), and Egypt (EG); Table S6: Parameters from OPLS-DA analysis for each comparison, using the ropls R Bioconductor package. For the analysis, the standard scaling method was used after log10 data transformation.